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Abstract. We study the motion of N — 2 overdamped Brownian particles in gravitational interaction in a 
space of dimension d — 2. This is equivalent to the simplified motion of two biological entities interacting 
via chemotaxis when time delay and degradation of the chemical are ignored. This problem also bears some 
similarities with the stochastic motion of two point vortices in viscous hydrodynamics [Agullo & Verga, 
Phys. Rev. E, 63, 056304 (2001)]. We analytically obtain the probability density of finding the particles at 
a distance r from each other at time t. We also determine the probability that the particles have coalesced 
and formed a Dirac peak at time t (i.e. the probability that the reduced particle has reached r = at 
time t). Finally, we investigate the variance of the distribution (r 2 ) and discuss the proper form of the 
virial theorem for this system. The reduced particle has a normal diffusion behavior for small times with a 
gravity-modified diffusion coefficient (r 2 ) = r 2 + (4kB /£/i)(T — T*)t, where feT, = Gmim2/2 is a critical 
temperature, and an anomalous diffusion for large times (r 2 ) oc t 1 ~ T * //T . As a by-product, our solution 
also describes the growth of the Dirac peak (condensate) that forms at large time in the post collapse 
regime of the Smoluchowski-Poisson system (or Keller-Segel model) for T < T c = GMm/(4fe). We find 
that the saturation of the mass of the condensate to the total mass is algebraic in an infinite domain and 
exponential in a bounded domain. Finally, we provide the general form of the virial theorem for Brownian 
particles with power law interactions. 

PACS. 5.20.-y Classical statistical mechanics - 05. 45. -a Nonlinear dynamics and chaos - 05.20.Dd Kinetic 
theory - 64.60.De Statistical mechanics of model systems 



1 Introduction 



Systems with long-range interactions have recently been 
the object of considerable interest [T]. One usually consid- 
ers isolated systems in which the particles evolve accord- 
ing to deterministic Hamiltonian equations. These sys- 
tems are described by the microcanonical ensemble. Ex- 
amples of such systems include self-gravitating systems, 
two-dimensional point vortices, the Hamiltonian mean 
field (HMF) model etc. However, one may also consider 
dissipative systems in which the particles, in contact with 
a thermal bath, evolve according to stochastic Langevin 
equations. These systems are described by the canoni- 
cal ensemble. The statistical mechanics of Hamiltonian 
and Brownian systems with long-range interactions is dis- 
cussed in P] at a general level. In this paper, we consider 
the case of Brownian particles in gravitational interac- 
tion. It is known that this system bears deep analogies 
with simple models of bacterial populations experiencing 



chemotaxis in biologjQ (see, e.g., Ref. [5] for a descrip- 
tion of this analogy). In a proper thermodynamic limit 
N — > +oo, the mean field approximation becomes ex- 
act and the dynamics of these systems is described by 
the Smoluchowski-Poisson system (gravity) [5] or by the 
Keller-Segel model (chemotaxis) [7|. These equations dis- 
play rich phenomena such as collapse and evaporation. In 
particular, in d = 2 dimensions, the evolution leads to 
the formation of Dirac peaks if the temperature is below 
the critical value ksT c = GMra/A [8|9I10] . These systems 
have been studied essentially in the mean field approxi- 
mation, i.e. for a large number of particles. The case of a 
finite number of particles can be studied numerically by 



These systems are isomorphic up to a change of notations. 
In this paper, we shall use the notations of astrophysics because 
they are closer to the notations that are familiar in physics 
and thermodynamics. However, our results can be transposed 
easily to the biological context. We refer to Perthame [3] for 
a complete bibliography of the chemotactic problem from the 
viewpoint of applied mathematics and to Chavanis & Sire for 
additional references in physics [4]. 
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solving the A-body stochastic equations. Numerical re- 
sults will be presented in a companion paper In the 
present paper, we consider the extreme case of only A = 2 
Brownian particles in gravitational interaction that can be 
solved analytically. We analytically obtain the probability 
density of finding the particles at a distance r from each 
other at time t and determine the probability that the 
particles have coalesced and formed a Dirac peak at time 
t. We also investigate the variance of the distribution (r 2 ) 
and discuss the proper form of the virial theorem for this 
system. In particular, we show that the virial theorem ob- 
tained in [T2] is only valid as long as the particles have 
not formed Dirac peaks. 

The paper is organized as follows. In Sec. [5] we recall 
the A-body coupled stochastic equations describing the 
evolution of self-gravitating Brownian particles and specif- 
ically consider the case A = 2. We introduce the center of 
mass and the reduced particle. We show that the center 
of mass undergoes a pure Brownian motion and that the 
reduced particle undergoes a Brownian motion in a cen- 
tral potential U = Gm\mi lnr. We also recall the "naive" 
virial theorem obtained in [12] and discuss, with a new 
light, the distinction between the critical temperatures 
fcsT c = Gmirri2/4 and ksT* = G?nim2/2. In Sec. [3] we 
study the motion of a Brownian particle (reduced parti- 
cle) in an attractive central potential U = GTO1TO2 hir in 
d = 2. We show that the corresponding Fokker-Planck 
equation is equivalent to a Schrodinger equation (with 
imaginary time) with a potential V = —D(a/r) 2 . This 
equation can be solved analytically in terms of Bessel func- 
tions. Then, we can obtain various analytical results such 
as the probability to find the reduced particle at position 
r at time t, the probability that the particle reaches the 
origin for the first time between t and t + dt, the prob- 
ability that the particle has reached the origin at time t 
and the variance (r 2 ) of the distribution. We find that 
the reduced particle has a normal diffusion behavior for 
small times with a gravity-modified diffusion coefficient 



(r 2 ) 



(4&b/£/z)(T — T*)t and an anomalous dif- 



fusion for large times 



x 



t 1 T '/ T . In particular, the 



variance increases with time when T > and tends to 
zero for t — > +oo when T < T„. In Sec. 2] we consider 
the case of two self-gravitating Brownian particles in a 
bounded domain and discuss the differences with the case 
of an infinite domain. Finally, in Sec. [5] we show that 
our study also describes the large time asymptotics of 
the Smoluchowski-Poisson system (or Kcller-Segel model) 
for T < T c = GMm/ (4k b)- Indeed, in the post-collapse 
regime, the system is made of a growing central Dirac peak 
(condensate) surrounded by a dilute halo whose dynam- 
ical evolution is eventually described by a Fokker-Planck 
equation similar to the one studied in the case of A = 2 
particles. We find that the saturation of the mass of the 
condensate to the total mass is algebraic in an infinite do- 
main and exponential in a bounded domain and we char- 
acterize it precisely. In Sec. [S] we briefly generalize our 
results to the logarithmic Fokker-Planck equation in d di- 
mensions. The Appendices provide complements such as 
the deterministic limit T = (Sec. 10]). the van Kampcn 



classification (Sec. 10, the correlation functions (Sec. [Gj 
and the general form of the virial theorem for Brownian 
particles with power law interaction (Sec. [XJ) . 

We may note that our study bears some similarities 
with the stochastic motion (induced by viscosity) of two 
point vortices studied by Agullo & Verga [13] ■ However, 
there also exists crucial differences between the two prob- 
lems since in our case the interaction is radial leading to 
the formation of Dirac peaks while in the case of point 
vortices the interactions is rotational leading to the for- 
mation of a spiral structure. 

We may also note that the statistical mechanics of 
A = 2 particles in gravitational interaction has been con- 
sidered by Padmanabhan [T3] in d = 3 (and generalized 
by Chavanis [H] for the dimensions d = 1 and d = 2) 
in the microcanonical and canonical ensembles. However, 
these authors consider the equilibrium statistical mechan- 
ics of A = 2 sclf-gravitating particles in a box, and with 
a small-scale cut-off, while we consider here the dynami- 
cal evolution of A = 2 self-gravitating Brownian particles 
in a finite or infinite domain without small-scale cut-off. 
Therefore, we address the time dependent problem and 
investigate the formation of Dirac peaks. 

Finally, the particular character of the dimension 
d = 2 in gravity is well-known. We refer for example to 
[15I16I14I17I18I19I20I9I12I2T] for more details and further 
references. 



2 The position of the problem 

2.1 The A-body problem 

We consider a system of A overdamped Brownian parti- 
cles with mass m a in gravitational interaction in a space 
of dimension d. Their motion is described by the coupled 
stochastic equations [12] : 

fly 1 

-£■ = —» Vatf(r 1 ,...,r JV ) + v^B a (t), (1) 

dt c,m a 



with 



C/(ri,...,r7v) 



G 



E 



d- 2 ^ |r Q -r g \ d - 2 '' 



(2) 



for d 2 and 

U(rx, ...,r N ) = G^2 m « TO /3 In l r a - r fs\, (3) 



a</3 



for d = 2. Here, £ is the friction coefficient and B Q (t) is a 
white noise satisfying (B Q (i)) = and (B a:i (t)Bpj(t')) = 
5ijS a p8(t — t') where a = 1, N refers to the particles 
and i = 1, d to the coordinates of space. The diffusion 
coefficient is given by the Einstein relation 



k B T 



(4) 
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where T is the temperature. We assume that the friction 
£ is the same for all the particles. 

From these stochastic equations, it is possible to derive 
the Fokker-Planck equation for the iV-body distribution 
Pjv( r i) —5 r JV) t) and then write the BBGKY-hierarchy 
for the reduced distributions |2I12) . Let us consider for 
brevity the single-species system. The proper thermody- 
namic limit corresponds to N — > +00 in such a way that 
the normalized temperature r\ = j3GMm/ R d ~ 2 is of order 
unity. In that limit, it can be shown that the mean field 
approximation becomes exact so that the iV-body distri- 
bution factorizes in a product of N one-body distributions: 
P/v (ri , r jv, £) = I1q -Pi ( r a ,t) [2j . Furthermore, the one- 
body distribution, or equivalently the smooth density field 
p(r,t) = NmP\(r,i), is solution of the Smoluchowski- 
Poisson system [2J: 



dp 
dt 



= V ■ 



1 

I 



k B T. 



(5) 



S d G P . (6) 

The equations generalizing Eqs. ©-(JB]) for the multi- 
species case are given in [12122] . 

Up to a change of notations, these equations are iso- 
morphic to a simplified version of the Keller-Segel model 
of chemotaxis that is valid in the limit of large diffusivity 
of the chemical and in the absence of degradation [5] . 

2.2 The case N ~ 2: the reduced particle 



Therefore, the center of mass undergoes a pure Brownian 
motion of the form 



with a diffusion coefficient 



D» = 



k B T 



(13) 



(14) 



Concerning the motion of the reduced particle, we have 
dr G . r 



= V2D 2 B 2 (t) - V2DiB!(i) = S(f), (15) 
where the noise satisfies 

(S i (t)S j (t')) = ^5 ij 5(t-t'). (16) 

Therefore, the reduced particle undergoes a Brownian mo- 
tion in a central potential of the form 



dr 



Gmim2 r 



with a diffusion coefficient 

k B T 



D 



(17) 



(18) 



From now on, we consider only N — 2 self-gravitating 
Brownian particles in d = 2. In that case, the stochastic 
equations ([I}- (J3j> reduce to 

dri Gm 2 ri - r 2 



dt i |r!-r 2 | 2 

dr 2 Gm\ Y\ — r 2 
~dt = £ |r-! — r 2 | 2 " 



h V2Di~Bi(t), (7) 
V^B^f), (8) 



with D\ = ksT/^mi and D 2 = kol '/ 't;m 2 . Like for the 
standard two-body problem, we introduce the center of 
mass 

roiri + m 2 r 2 
R= — , M = mi+m 2 , (9) 

and the reduced particle 



r = r 2 - n, 



m\m 2 



mi + m 2 

Concerning the motion of the center of mass, we have 
rfR 1 



(10) 



dt M 
where the noise satisfies 



(m lv /2AB 1 (i)+m 2V /2Z^B 2 (i)) = Q(f), 



(11) 



(Qi(t)Qoit')) = ^S ij 6(t-t'). (12) 



2.3 The naive virial theorem 

Let us introduce the total moment of inertia 

Itot(t) = m 1 (r 2 1 )+m 2 (r 2 2 ). (19) 

In [T2] (see also Appendix an exact closed expression 
of the virial theorem valid for an arbitrary number of self- 
gravitating Brownian particles in d = 2 has been obtained. 
For N = 2, it writes 



^Itot = 2fc B (T-T c ), 



with the critical temperature 



kuT r 



Gm\m 2 



(20) 



(21) 



It is instructive to recover this result in a different manner. 
The positions of the particles 1 and 2 can be expressed in 
terms of r (reduced particle) and R (center of mass) as 

MR-m 2 r MR + TOir , , 
r i = ) r 2 = — • ( 22 ) 



M 



M 



Substituting these relations in Eq. (fT9"]) . we obtain after 
straightforward algebra 



Itot(t) = M(R 2 )+p(r 2 }, 



(23) 
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a relation which was of course expected. Now, the Fokker- 
Planck equation associated with the stochastic motion 
(fTT|) of the reduced particle is 



? at 



knT Gmim? r , 
— VP + P — . 24 



Taking the time derivative of 



= / Pr 2 dr, 



(25) 

and using simple integrations by parts, we naiveljQ obtain 

(26) 



This relation exhibits a critical temperature 

Gmim 2 
krT* = . 



(27) 



2.4 The problem 

In fact, there is a flaw in the above derivation of the virial 
theorem because we have naively assumed that the nor- 
malization J P(r, t) dr = 1 is conserved in time. However, 
as we shall see, this is not correct. The normalization is 
not conserved in time because the reduced particle can 
reach the origin r — and be "lost" by the system (if it 
reaches the origin, it remains there for ever). This corre- 
sponds to the coalescence of the two particles, resulting in 
the formation of a Dirac peak, i.e. a new particle of mass 
mi + n%2- As a result of these "trapping" events 

Jp{v,t)dv^l, (33) 

and we must reconsider the problem in more detail. 



3 Brownian particle in a Newtonian potential 
in two dimensions 



Introducing the moment of inertia of the reduced particle 
I(t) = n(r 2 ), we can rewrite Eq. (|26p as 



(28) 



The mean square displacement of the reduced particle sat- 
isfies 

(r 2 ) = 4fcB 2 



This is a normal diffusion with a gravity modified diffusion 
coefficient 

0(11 -l£ ('-£)• (30) 

The variance increases for T > and tends to zero in 
a finite time for T < T*. On the other hand, the Fokker- 
Planck equation associated to the stochastic motion (|13l) 
of the center of mass is simply 



? at m 

and we classically obtain the relation 



V dt 



k R T. 



(31) 



(32) 



Finally, summing Eqs. (|26p and (|32p and using Eq. (|2l}|) . 
we recover Eq. (|20[) . We now clearly see the origin of the 
two temperatures T c and = 2T C that were reported in 
[T2"] . In the case N = 2, the critical temperature T* is as- 
sociated to the dynamics of the reduced particle while the 
critical temperature T c enters in the expression of virial 
theorem for the total moment of inertia (reduced particle 
and center of mass). This distinction is further discussed 
in Appendix lAl in the general case of TV particles. 



We shall see later that this expression is in fact incorrect. 



3.1 The Fokker-Planck equation 

Let P(r, t) denote the probability density of finding the 
reduced particle in r at time t. The evolution of P{r,t) is 
governed by the Fokker-Planck equation^ 



? at 



knT Gm^m-y r , 

VP + P — — |. (34) 



The initial distribution Po( r ) is normalized such that 
/ Po( r )^ r = 1- Introducing 



U = Gm,\m2 lnr, 



/i£ k B T 
the Fokker-Planck equation can be rewritten 
dP 

— = V- [D (VP + PPVU)]. 



(35) 



(36) 



In the absence of small and large scale cut-offs, this 
equation has no steady state since the distribution P = 
A/r° mim2 is not normalizablc. We assume that the 
initial distribution Pq(t) is radially symmetric, so that 
P(r, t) is radially symmetric for all times. Therefore, we 
can write the Fokker-Planck equation in the form 



dP _ 1 d 
dt r dr 



M 



dP BGmimo 
— + P 1 - — 

or r 



(37) 



As discussed previously, the probability is not conserved 
because the reduced particle may reach the origin r = 
and form a Dirac peak (the two particles coalesce). The 
probability that the particle has not reached r = at time 
t is 



r+oo 

X(t) = / P{r,t)2iTrdr. 
Jo 



(38) 



Note that, for d = 1, the Fokker-Planck equation for the 
reduced particle corresponds to a V-shaped potential U(x) = 
Gmim 2 \x\ that relaxes to P a [x) = ^Gmim 2 e fGmim2111 [53]. 
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Taking the time derivative of this quantity and using the 
Fokker-Planck equation ([37]) we obtain 



(39) 



X(t) = -2<KDj3Gm 1 m 2 P(0,t), 



which is non zero since P(0,t) 7^ 0. Therefore, the proba- 
bility for the particle to form a Dirac peak between t and 
t + dt (i.e. to reach r = for the first time between t and 
t + dt) is 

Xo{t) = 2nDpGm 1 m 2 P(0, t), (40) 

and the probability for the particle to have formed a Dirac 
peak at time t (i.e. to have reached r = at time t) is 



X D (t) = 2nD(3Gm 1 m 2 [ P(0,t)cLt. 

Jo 



(41) 



We obviously have Xd($) = 1 — x(t)- We can now obtain 
the proper form of the virial theorem associated to the 
Fokker-Planck equation ([34|) . Introducing the moment of 
inertia of the reduced particle 



I{t) = J P^ir 2 dr, 



we easily obtain the virial theorem 



(42) 



(43) 



instead of Eq. (|28[1 . It has to be noted that this relation is 
not closed since it depends on x{t) that must be obtained 
by solving the Fokker-Planck equation (|34)) . 



For a spherically symmetric distribution, the Schrodingcr 
equation (|47p can be rewritten 



with 
V(r) 



1 a 1 d ( d-iTM \ dU 



Making the separation of variables 
V = e" A V(r), 
we obtain the eigenvalue equation 
1 d 



(50) 
(51) 



r d 1 ^ r 



r d - l D(r) 



dr 



+ (V{r) + \)0 = O. (52) 



Let us note X n the eigenvalues and <f> n the corresponding 
eigcnfunctions. The eigenfunctions are orthogonal with re- 
spect to the scalar product 



(f\9) = / /(r) 3 (r)5 rf r d " 1 dr. 



We also normalize them so that (4> n \4> m ) = S n 
any function can be expanded in the form 



+00 



h(r) = J2(h\MMr 



(53) 
Then, 

(54) 



3.2 The associated Schradinger equation 

Let us consider a general Fokker-Planck equation of the 
form 



8P_d_ 
dt dr 



D[ -r r +PP ir 
\ or or 



(44) 



For a spherically symmetric distribution in d dimensions, 
it can be rewritten 



dP 
~dt 



1 



(45) 



As is well-known [53] , we can transform this Fokker-Planck 
equation into a Schrodingcr equation (with imaginary 
time) by setting 



P(r,t) 



e 2 



V>(r,t). 



(46) 



This yields 



dtp d ( dip 
dt dr \ dr 

with the potential 



1 „ d 



V{r)tP = (H + V)ip, (47) 



dU 



^=2^'l^J-4^ 2 (^ 



dU 



(48) 



If the spectrum is continuous, the sum over n must be 
replaced by an integral over A > 0. 



3.3 The general solution 

We consider the Green function P(r,t\ro) which corre- 
sponds to the initial condition 



(55) 



The solution on the Fokker-Planck equation (|4"5|) can be 
expanded on the eigcnfunctions in the form 



+00 



P(r,t\r ) = ^2A n e- Xnt e-^ u ^Un(r). (56) 

n=l 

Noting that 



„ -j^T = H <Pn(ro)(t>n(r), 
and using the initial condition (|55p . we finally obtain 



+00 



P(r,t\r )=e 



- p -kP(U{r)-U{r )) ST P -A„t 



J2 e ~ Xnt Mro)K(r). (58) 
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3.4 The case of a logarithmic potential in d = 2 

For the Fokkcr-Planck equation ([57]) . we have d = 2, 
D = k B T/(£_fi) and U = Gm 1 m 2 \nr. The potential V(r) 
arising in the corresponding Schrddinger equation (1491) is 



V(r) = -D 



( j3Gmim2 



\ 2r 

Therefore, if we assume that initially 

S(r - r ) 



P(r,0) 



27rr 



(59) 



(60) 



the solution of the Fokker-Planck equation (|37|) can be 
written 



P 



(»".*) = (7) I e- A VA(r o )0A(r)dA, (61) 



where <fi\(r) is solution of the differential equation 

A 



where 



rV + rtfif + [ —r 2 - a 2 ) <f> = 0, 



/3Gmim 2 

2 ~ T 7 ' 



(62) 



(63) 



Equation (|62|) is a Bcssel differential equation that can 
be solved analytically. The solutions that are finite at the 
origin are of the form 



Mr) = Ja(VX/Dr). 



(64) 



Substituting Eq. (l64l) in Eq. (juTj) . the solution of the 
Fokkcr-Planck equation (|3"T)) can be written 

P(r, i) = 2£ (^) ° A 2 y + °° e- m2 * J a (\r ) J a (Ar) ArfA, 

(65) 

where we have made the change of notation A — > DX 2 for 
convenience. The normalization constant A is determined 
so as to recover the initial condition (|60|) as t — > 0. Taking 
t = in Eq. ([6"5]) . wc get 



P{r,0)=2D(^Y A 2 £ J a (Xr )J a (Xr) XdX. 



(66) 



Using the closure relation 

»+oo 



J a (ux)J a (vx)x dx = —5(u — v), (67) 



wc obtain 



on 42 

P(r,0) = <f(r-r„). 

TO 



(68) 




T=T,/2 



Fig. 1. Probability density of finding the particles at a dis- 
tance r from each other, at different times t and for a tem- 
perature T = T*/2 (a = 2). From top to bottom: t — 
0.01,0.025,0.05,0.1,0.2,1. 



Comparing with Eq. $U$, we find that A 2 = 1/(4ttD). 
Therefore, the solution of the Fokker-Planck equation (f3"T|) 
with the initial condition (l60l) is 



1 /r \ a 



+00 



e- DA2t J a (Ar )J a (Ar) XdX. 



(69) 



Using the identity 

r*+co 







_ 2 2 1 ° 2 +f 2 /a/3 

e p J 1 (ax)J 1 ((3x) xdx = — -e 4 ? 2 I 7 — r 

2p^ \2p 



(70) 



valid for 7 > — 1, we find that it can finally be written 

p M ) - (2) 



e — J n — M. (71) 

47rL>7i a \2Dt) y ' 



The distribution P(r, i) is plotted in Fig. [T] at different 
times and for T/T* = 1/2 (corresponding to a = 2). 
Using the identity 



Ia(x) ~ 

we get for t —> 0: 

P(r,i) 



(x — > +00) 



(72) 



(73) 



which tends to Eq. (|6"0"|) as expected. On the other hand, 
for t — > +00, the probability tends to zero meaning that 
the particle has been absorbed in r = after a sufficiently 
long time so that it is ultimately lost by the system. For 
r — > +00, using the identity (|72[). we have 



p ( M) r°y 1 - 1 

V r ) 47r y'TirroDt 



For r = 0, using the identity 

/q(x) ^ P(o+1) (2 



(a: -> 0) 



(74) 



(75) 
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we get 



K 1 ' \2J 47rr(a+l) {Dt) l + a 



e 4Dt 



(76) 



Finally, for ,8 = 0, the probability density (|TTj) becomes 
1 



P(r,f) 



2 , 2 
+ 



— — e 4Dt in 

AnDt \2DtJ 



(77) 



which is the solution of the diffusion equation in d = 2. 
Indeed, in the limit of infinite temperature, the gravity is 
negligible with respect to diffusion. 



3.5 The probability to form a Dirac peak 

The probability that the particle has not reached r = at 
time t is given by Eq. ([55)1 . For the distribution (|TTj) . the 
integral can be performed analytically and we obtain 



where 



x(t) 



r a (x) 



4Dt 



r(a) 



t" 



"'dt, 



(78) 



(79) 



is the incomplete Gamma function. The probability decays 
because, as time goes on, the particle has more and more 
chance to reach r = and form a Dirac. The probability 
that the particle reaches r = for the first time between 
t and t + dt is given by Eq. (|40|) . Combining this relation 
with Eq. ([To]) , we obtain 



XD 



D 



r \ 2a 1 



1 



r(a) (Dt) 



l + Q 



e ro* . 



(80) 



Integrating Eq. (|80D . we obtain the probability that the 
particle has formed a Dirac peak at time t: 



Xn(t) = 



r, 



a \4Dt J 

7>T 



(81) 



and we check that x(t) = 1 — Xo(t) as expected. The 
evolution of \d {t) for different values of the temperature 
is shown in Figs. [2] and |3] 

For t — ¥ 0, using the expansion 



r a {x) 



we obtain 

XD(t) > 



1 



r(a) \ADt J 



(x^+oo), (82) 



e _ 35t, (t->0). (83) 



We see that the probability x_o(t) — > for t — >• due to 
the exponential factor. This tendency is reinforced by the 
algebraic factor for T > T* (a < 1) while it is reduced for 
T <T„ (a> 1). 




Fig. 2. Evolution of the function xo{t), giving the probability 
that the particle has formed a Dirac peak at time t, for different 
values of the temperature (we have taken T/T» = 1/2, T/T* = 
1 and T/T» = 2). The probability increases more rapidly at 
smaller temperatures. 




Fig. 3. Same as Fig. [2] for larger times. 



For t — > +oo, using the expansion 
.r" 



we obtain 

XD® = 1 



(x -> 0), 



(84) 



or(a) \4i)i 



(i^+oo). (85) 



Therefore, the probability that the particle has not formed 
a Dirac at time t decreases algebraically as t~ a . Equation 
13 can be written in the form 



where the time 



XD(t) = 1 



t*(a) 



t*(a)\ a 



AD 



[ar(a) 



l/a ' 



(86) 



(87) 



gives an idea of the rapidity at which the Dirac forms as 
a function of the temperature a = T*/T. The function 
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Fig. 4. Evolution of t*(a) as a function of the temperature 
a = T»/T. 




Fig. 5. Evolution of as a function of the temperature 

a = T»/T. 



i* (a) is represented in Fig. S] and its asymptotic behav- 
iors are given in Appendix [Ej We find that 4_Di* (a) /r 2 ~ 
1.78107... for a = and 4Dt»(o)/rg ~ 2.7182818/a... for 
a — > +oo. We note that i*(a) tends to a finite value for 
a — > while we know that the system does not form a 
Dirac peak for a = (indeed 1 — (t*(a)/t) a — s- for a — s- 0). 
Therefore, the physical interpretation of t*(a) should be 
considered with care. Another measure of the effect of the 
temperature on the formation of the Dirac is provided by 
the quantity 



Xd(1) = 



no)' 



corresponding to X£>M evaluated at t = r^/AD. This func- 
tion is represented in Fig. [5] and its asymptotic behaviors 
are given in Appendix|El We find that X-d(I) ~ 0.219384a 
for a — s- and 1 — X-d(I) ~ a a+i/2 f° r a +°°- 

Finally, the normalized probability density can be writ- 
ten 

P tot (r,t) = P(v,t)+ X D(t)6(r), (89) 



where P(r,t) is given by Eq. §TB) and xd(*) by Eq. ([STjl. 
We readily check that / P*ot(r, t) dr = 1. 



3.6 The moment of inertia 

The moment of inertia of the reduced particle is defined 

by 

r+oo 

I(t)= / P(r,t)fir 2 2nrdr, (90) 

and the variance of the distribution (mean square displace- 
ment) is 

<r 2 ) = ^. (91) 

For the density distribution given by Eq. (|71l) , the integral 
can be calculated explicitly yielding 

I(t ) = JH-(A-Y Dt e-& 
(t) r(a) \4Dt ) 



+ti(rl+4D(l-a)t) 



A.) 



r(a) 



(92) 



In Appendix [B] we check that this relation is consistent 
with the virial theorem (|4"3"|) . The evolution of the moment 
of inertia is represented in Fig. [5] for different values of the 
temperature. For T > T*, the moment of inertia increases, 
for T = T* the moment of inertia is constant I(t) = \ir\ 
and for T < T* the moment of inertia decreases. This will 
become clear from the asymptotic behaviors. 
For t —> 0, we get 



I{t) ~ nrl+4Dn{l-a)t, 



which can be rewritten 



CM 



(93) 



(94) 



In that case, we have a normal diffusion with a gravity- 
modified diffusion coefficient 



D(T) = -£.(T-T,). 



(95) 



For T > the variance increases with time while for T < 
T* it decreases. This expression agrees with the naive virial 
theorem (f2"9"| . Indeed, for small times, the probability for 
the particle to reach r = is exponentially small so that 
the probability is conserved and the naive virial theorem 
holds since there is no Dirac peak. 

For t —> +oo, using the expansion (|84p. we get 



This corresponds to 



ar(a) 



(96) 



(97) 



For T > T*, i.e. a < 1, the variance increases and goes to 
+oo for large times. In that case, we have an anomalous 
diffusion (r 2 ) ~ t a with an exponent a = 1 — T*/T. The 
evolution is always sub-diffusive. The origin of the anoma- 
lous diffusion is related to the fact that the particle can 
be trapped at r = (and form a Dirac peak). For T < T„, 
the variance decreases and goes to for t — > +oo. 
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T>T, ^^^^ 




T = T, 


T<T, 

i ; i ; 1 1 1 , 



10 



4Dt/r„ 




Fig. 6. Time evolution of the moment of inertia for different Fig. 7. Time evolution of rp for different values of the temper- 
values of the temperature (we have taken T/T, = 1/2, T/T, = ature (we have taken T/T, = 1/2, T/T, = 1 and T/T, = 2). 
1 and T/T, = 2). 



3.7 The most probable position 



The most probable value rp(t) of the distribution P(r, t) is 
obtained by maximizing P(r,t), or cquivalently \nP(r,t), 
with respect to r. This gives 

2aDt , r P I' a (^) 



rpr r I a (^) 
Using the recurrence relation 



r a {x) = i a+1 {x) + -i a (x), 



we obtain 



7"0 



\ 2Dt 
J I rprg \ 
±a \ 2DI I 



(99) 
(100) 




This equation can be rewritten in the parametric form 
r P r 2Dt I a +i(x) 



~ 2Dt ' rl xl a (x) ' 

or equivalcntly 

r P _ I a +i{x) 2Dt _ I a +i(x) 



Ia{x) 



xl a (x) ' 



(101) 



(102) 



Fig. 8. Time evolution of r» for different values of the tem- 
perature (we have taken T/T, = 1/2, T/T, = 1, T/T, = 2 and 
T/T, = 4). 



Therefore, the radius rp(t) is decreasing for any temper- 
ature and it goes to zero in a finite time t c {a) depending 
on the temperature. Some curves are represented in Fig. 

m 

The most probable value r, (t) of the radial distribution 
P{r,t)r is solution of 



2(a- l)Dt 



jl I r,r a 
a \ 2Dt , 



which gives rp(t). Using the asymptotic expansions of 
I a (x), we find that 



r,r Q 

Using Eq. (|99]l. we get 



(106) 



and 



ro 



r JL = 1 

ro 



a + 2 
a+l 



2a + 1 2Dt 



1/2 



(t -»> o), 



where 



Vi - t/t c , (t -> t c ), 



2Dt c _ 1 
~rf~ 2(a + l)- 



(103) 

(104) 
(105) 



r. 



2Dt I a+1 (fgf) 



r*r /o(i^) 
This equation can be rewritten in the parametric form 
r,r 2Dt _ I a +i(x) 1 



(107) 



x — 

2Dt' 

or equivalcntly 

r* _ I a +i(x) 



Ia(x) 



2Dt 



xl a (x) x- 

= 4+1 (x) 
a;/ (x) 



(108) 



4' ( 109 ) 
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which gives r*(t). Using the asymptotic expansions of 
I a (x), we find that 



±=1 

ro 



2a -1 2Dt 



(t-^0), 



and 



(t -> +oo). 



(110) 



(111) 



Therefore, the radius (t) is always increasing for T > 
2T*. For T < 2T*, it starts to decrease before finally in- 
creasing. For T = 2T*, we can use the identities 



h/2(x) 



h/2(x) = \ — 



cosh(a-) 



sinh(a;), 



sinh(x) 



(112) 



(113) 



so that r*/ro = 1/ tanh(x)— ljx = L(x) (where L(x) is the 
Langevin function) and 2Dtjr\ — l/(a;tanh(x)) — l/x 2 . 
Using the asymptotic expansion of tanh(x) for x — > and 
x —> +oo, we find that 



l + 2e-Bt, ((^0). (114) 
Some curves are represented in Fig. [8] 

3.8 Reflecting boundary conditions 

In the previous sections, we have considered the case of 
absorbing boundary conditions at r — 0. They lead to the 
formation of a Dirac peak with an amplitude xo(t) grow- 
ing with time. We shall now consider the case of reflecting 
boundary conditions and determine their domain of ex- 
istence. Reconsidering the calculations of Sec. 13.41 these 
boundary conditions correspond to solutions of the form 



<f>x(r) = J- a (VX/Dr), 



(115) 



that diverge at the origin. Substituting Eq. f|l 15|) in Eq. 
(jUTj) and repeating the calculations of Sec. 13.41 we obtain 

According to identity (|70p . this solution is valid only for 
a < 1 (i.e. T > T*). This is confirmed by considering the 
equivalent of P(r, t) close to r = 0: 



P(r,t) 



1 



4" 



1 



r f3Gm im2 4 n r(l - a) (Dt) l - a 



(117) 



This distribution diverges at the origin and it is normaliz- 
ablc iff T > T*. Since the distribution P(r,t) diverges at 
the origin, we must replace Eq. (|39p by 



xW = -2^hm,f^ + P^™ 
r-s-o V or r 



(118) 



T=2 T, 




Fig. 9. Probability density of finding the particles at a dis- 
tance r from each other, at different times t and for a tem- 
perature T = 2T» (a = 1/2). From top to bottom: t = 
0.01,0.025,0.05,0.1,0.2,1. 



Using Eq. f|116[) and the expansion 



Ia{x) 



r(a+l) 



x\ a fx\ a + 2 

2) +7^(2) + -' ( 119 ) 



r(a + 2) 



valid for x — > 0, we find that \ = 0. Therefore, the nor- 
malization is conserved in time and there is no Dirac peak 
formation. 

These results arc consistent with the van Kampen clas- 
sification of singularities (see Appendix [F]) . For T < 
(i.e. a > 1 or f3Gmim2 > 2), the singularity at r = 
behaves as an adhesive boundary. In that case, the solu- 
tion is unique and no boundary condition has to be fixed 
by hand. It is given by Eq. (fTTj) leading to a Dirac peak 
(x ^ 0). On the other hand, for T > (i.e. a < 1 or 
(3Gmini2 < 2), the singularity at r = behaves as a regu- 
lar boundary. In that case, the boundary condition can be 
absorbing, reflecting or mixed. It has to be fixed by hand. 
In the previous sections, wc considered a purely absorbing 
boundary condition and, in the present section, we con- 
sidered a purely reflecting boundary condition. In fact, for 
a < 1 (i.e. T > T„), the general solution can be written as 
a "mixture" of the two previous solutions 



P(r,t)=vP+(r,t) + (l-v)P-(r,t), 



(120) 



where P+(r,t) is the distribution (fTTj) . P-(r,t) is the dis- 
tribution (|116[) . and v is a parameter taking values be- 
tween v = 1 (purely absorbing) and v = (purely reflect- 
ing). 

Let us specifically consider the distribution (| 1 16[) . It 
is plotted in Fig. [9] at different times for T/T* = 2 (cor- 
responding to a = 1/2). Substituting Eq. (|116[) in the 
expression (|90[) defining the moment of inertia of the re- 
duced particle, and carrying out the integrations, wc find 
that 

I(t) = ^rg + 4D/x(l - o)t. (121) 



This expression agrees with the naive virial theorem 
Indeed, in the present case %(i) = 1 since there is no Dirac 
peak formation, and the exact expression (|4"3"j) of the virial 
theorem reduces to Eq. 
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2Dt/r 2 

Fig. 10. Time evolution of r P for T = 2T» (a = 1/2). The 
upper branch corresponds to the position of the local maximum 
and the lower branch to the position of the minimum. 




Fig. 11. Time evolution of r» for different values of T (corre- 
sponding to a — 0.3, 0.5 and 0.7). In the case 1/2 < a < 1, the 
upper branch corresponds to the position of the local maximum 
and the lower branch to the position of the minimum. 



On the other hand, the most probable value rp(t) of 
the distribution P(r,t) is obtained by maximizing P(r,t), 
or equivalenly In P(r, t) with respect to r. In fact, since this 
distribution diverges at r = 0, the global maximum (infi- 
nite) is rp(t) = 0. However, for sufficiently short times, the 
distribution P(r, t) also admits a local maximum r^ (t) 

and a local minimum r p (t). Proceeding as in Sec. 13.71 
they are the solutions of the equation 



AaDt r P 
r P r Q r 



J I TpTg 

ll ~ a \ 2Dt , 
j I rpro \ 
1 ~ a \ 2Dt ) 



(122) 



Setting x — -^m 
parametric form 

rp _ h-a(x) 
r a I- a {x) 



rpr ° this equation can be rewritten in the 



2a 2Dt _ h- a (x) 2a 

X ' 7-q xl- a (x) X 2 ' 



(123) 



which gives rp(t). Using the asymptotic expansion of I a (x) 
for x — > +oo, we find that 



' p 
ro 



= 1 - 



2a + 1 2Dt 



(t->0). 



(124) 



Setting x = this equation can be rewritten in the 

parametric form 



I-a(x) 



2a -I 2Dt _ h- a {x) 2a -I 

X Tq xl- a (x) X 2 



(126) 

which gives (t). Using the asymptotic expansion of I a {x) 
for x — > +oo, we find that 



1 - 2a 2Dt 



For 1/2 < a < 1, the radius r* 1 ^) of the local maximum 
is decreasing and it disappears at a time t en d(a) which 
depends on the temperature. For < a < 1/2, the radius 
r^\t) increases initially and behaves for large times as 



(f ^0). 



(127) 



r* - V 2 (l - 2a)Dt, (t -> +oo). (128) 
For a = 1/2, we can use the identities (|112[) and 



The radius rp (t) of the local maximum is decreasing for 
any temperature and it disappears at a time t enc i(a) which 
depends on the temperature. The curve corresponding to 
a = 1/2 is represented in Fig. [TOl 

Let us now consider the most probable value (t) of 
the radial distribution P(r,t)r. This distribution diverges 
at the origin r = for 1/2 < a < 1, vanishes at the origin 
for < a < 1/2 and is finite at the origin for a = 1/2. In 
the first case, the global maximum (infinite) is at rp = 
but, for sufficiently short times, the distribution P(r,t)r 

also admits a local maximum rp and a local minimum 

( 2 ) 

r p . Proceeding as in Sec. 13. 71 they are the solutions of 
the equation 



I-i/2(x) = \ — cosh(x), 



(129) 



so that r*/ro = tanh(x) and 2Dtjr\ = taiih(x)/x. Using 
the asymptotic expansions of tanh(x) for x — > and x — ) 
+oo, we find that 



— ~ 1 — 2e Dt 



(*->(>), 



1/2 



2D 



(130) 



(131) 



2(2o- l)Dt 



'l-o 



> 2Dt , 



r*r 



ro 



I-r 



, 2Dt , 



(125) 



establishing t en d = t c = ^pj for a = 1/2- Some curves are 
represented in Fig. 1111 
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4 The case of a bounded domain 

In the previous sections, we considered the case where the 
particles are free to move in an infinite domain. We shall 
now consider the case where the particles are confined in 
a bounded domain. As an idealization, we shall consider 
that the reduced particle evolves in a spherical box of ra- 
dius R. We first look for the existence of an equilibrium 
state. The steady solution of the Fokkcr-Planck equation 
satisfies 



VP e „ + /3G mi m 2 P t 



0. 



which is integrated into 



e 1 r l3Gm 1 m 2 



(132) 



(133) 



It is clear that this distribution is normalizable iff 
/3GmiTO2 < 2 so that an equilibrium state exists iff 
T > T*. For T < T», the singularity at r = is ad- 
hesive and a Dirac peak grows, ultimately absorbing the 
particle with probability one. For T > T*, the singularity 
at r = is regular and its nature (absorbing, reflecting 
or mixed) has to be fixed by hand. In the case of an ab- 
sorbing boundary, we ultimately obtain a Dirac peak of 
amplitude X-d(+°°) = 1, so that there is no equilibrium 
state outside the Dirac. In the case of a reflecting bound- 
ary, there is no Dirac peak and the system reaches an 
equilibrium state, given by Eq. (|133[) . with normalization 
/ P eq dr = 1. In the mixed case, we ultimately obtain a 
Dirac peak with amplitude %d(+oo) = v and an equi- 
librium state outside the Dirac, given by Eq. (|133[) . with 
normalization J P eq dr = 1 — v. 

Let us now solve the Fokker-Planck equation in 
a bounded domain. In that case, we need to impose the 
boundary condition 



or or 



0, 



(134) 



in r = R meaning that there is no flux of probability at 
the boundary. In terms of the eigcnfunction <f> defined in 
Eqs. (|4"6")l and (|5ip . this can be rewritten 



For the potential U = Gm\m 2 lnr, we obtain 
4>'{R) _ /3Gm im2 



2R 



(135) 



(136) 



This boundary condition implies that the eigenvalues are 
quantized. 

Let us first consider purely absorbing boundary condi- 
tions at r = 0. Making the change of notations A — > DX 2 , 
the eigenfunctions that are solution of Eq. (|6"2"j) and that 
are finite at the origin are given by 



<j) n (r) = A n J a (X n r). 



(137) 



Substituting this solution in Eq. (|136p . we find that the 
eigenvalues are determined by 



X n RJ' a (X n R) 

t- r — = —a. 

J a \X n R) 

Using the recurrence relation 

J' a (x) = J a -l(x) - ~Ja(x), 
X 

the foregoing equation can be rewritten 

Ja-l{X n R) = 0, 



(138) 



(139) 



(140) 



with X n 0. Therefore, the eigenvalues RX n (a) are the 
zeros of J a —i(x). Using the general identity 



J 2 (Xr)r dr 



R 2 



J' a {XRf 



1 - 



X 2 R 2 



J 2 a (XR) 



(141) 



together with the relation (|138[) . we obtain 



J J 2 (X n r)rdr = ^J 2 (X n R). (142) 



Therefore, the normalized eigenfunctions are 

1 



(t>n(r) = 



^RJ a {X n R) 



J a (X n r). 



(143) 



Finally, using Eq. (|5"8"|) . the solution of the Fokker-Planck 
equation (|34[) in a bounded domain can be written 



<t>n{ro)(f> n {r). 



(144) 



We can now redo the preceding analysis except that the 
results will be less explicit since they will be expressed in 
the form of series. 

The probability Xd($) for the particle to have formed 
a Dirac peak at time t is given by Eq. (|4ip . Using Eq. 
(I144[) and the equivalent 



J n (x) 
we obtain 



1 



r(n+l) \2 

P(0j t) = J_ ( r Ji) a —±—Yx 

y ' ' irR 2 V 2 J r(a+l) ^ 



(145) 



a 

JaiX 

n^O) —D\ 2 t 



J'a(^nR) 



(146) 



Therefore, the probability for the particle to have formed 
a Dirac peak at time t can be written 



X 



D (t) = l-^fl n 



— D\z t 



(147) 



with 



B n — 



2(3Gmim 2 fr \ a 1 xa -2 J a(X n r ) 



R 2 



r(a+l) K J 2 (X n R) 



(148) 
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For t 



-oo, wc obtain 



-DAT* 



(149) 



so that the probability that the particle has not formed a 
Dirac at time t decreases exponentially as e - m i* instead 
of algebraically in an unbounded domain (see Sec. I3.5[) . 
The exponential decay is controlled by the first eigen- 
value Ai(T) (fundamental) of the Schrodingcr equation. 
As we have seen, Xi(a)R is the first zero of J a _i(x). This 
is valid for any temperature. It is instructive, however, to 
determine the asymptotic behavior of Ai for a — > (i.e. 
T — > +oo). In that limit, Ai — > 0. Substituting the expan- 
sion 



J n {x) 



1 



1 



r(n+l) i> + 2) 



(150) 



valid for x — >• in Eq. (|140l) . we obtain 

2GTOim2 



(AiJ?)^ 



4a. 



On the other hand, in this limit, B\ 
(|149|) becomes 



XD (t)~l-e 
Using Eq. (|35[) . it can be rewritten 



2£Gmjm2 
k H TR 2 



Xn(t) 



1 



(151) 

1. Therefore, Eq. 

(152) 

(153) 



For T = (i.e. a = 1), we find that Aii? = jo,i where 
jo,i ^ 2.40482... is the first zero of J (x). Finally, for T -> 
0, we find that 



Aii? 



Gm\m,2 
2k B T ■ 



(154) 



Let us now consider purely reflecting boundary condi- 
tions at r = 0. Making the change of notations A — > DA 2 , 
the eigenfunctions that are solution of Eq. (|62j) and that 
diverge at the origin are given by 



(f> n (r) = A n J- a (\ n r). 



(155) 



Substituting this solution in Eq. (|136|) . we find that the 
eigenvalues are determined by 



X n RJ'_ a (X n R) 
J-a{X n R) 

Using the recurrence relation 

J'a( x ) = ~Ja+l(x) + ~Ja(x), 
X 

the foregoing equation can be rewritten 

Jl-a(KiR) = 0, 



(156) 



(157) 



(158) 



with A„ 7^ 0. Therefore, the eigenvalues RX n (a) are the ze- 
ros of Ji_ a (x). Using the general identity (|141[) together 
with the relation (|156[) . we find that the normalized eigen- 
functions are 

y/ltRJ- a (\ n R) 

This expression is valid for A„ 7^ 0. We must also add 
the eigenmode corresponding to Ao = whose normalized 
expression is 



Mr) 



1 - a 1 R 



TT R 



(160) 



Finally, using Eq. (|5"5)) . the solution of the Fokker-Planck 
equation (|34j) in a bounded domain can be written 



P(r,t) = P eq {r) + P) ^e-^Vn(r )^n(r), (161) 



where 



P eq (r) 



1 - a (R 



2d 



(162) 



ttR 2 \ r 

is the equilibrium distribution and the series run over n > 
1. Using Eq. (|161|) and the equivalent (|145l) . wc obtain for 

r ->• 0: 



P(r,i)-P e9 (r) 



J 1 (2r ) Q 

Tri? 2 r 2a r(l - a) 



The distribution diverges at the origin and it is normaliz- 
ablc iff a < 1, i.e. T > T*. On the other hand, using Eqs. 
(|16ip and pi8p . we find that x = so that there is no 
Dirac peak in that case. For t — > +00, we get 

P(r,t) - P eq (r) ~ ^y M r )Mr)e- DX ^, (164) 

so that the distribution converges exponentially rapidly 
towards the equilibrium state as e~ DXlt . The exponential 
relaxation time is controlled by the first eigenvalue Ai (a) 
(fundamental) of the Schrodingcr equation. As we have 
seen, i?Ai(a) is the first zero of J\- a {x). This is valid for 
any temperature T > T„. For T — > +00 (i.e. a — > 0), 
we find that Aii? — > j\ t \ where ~ 3.83171... is the 
first zero of J\{x). For T = (i.e. a = 1), wc find that 
Aii? = jo,i where 70,1 — 2.40482... is the first zero of Jo(x). 

Remark: wc could also consider the case of N = 2 
self-gravitating Brownian particles with a short-range reg- 
ularization. This could be due to a softened potential 
U = Gmi7B2 ln(\/r 2 + e 2 ), to a hard core a or to an ex- 
clusion principle such as the Pauli exclusion principle for 
fermions in quantum mechanics. In that case, there is no 
Dirac peak formation. There exists an equilibrium state 
for any temperature in a box and for T < in an infinite 
domain j!2j . For low temperatures, the two particles are 
at a typical distance a from each other so that the equilib- 
rium state is controlled by the small-scale cut-off. For high 
temperatures, the particles are at a typical distance i? (the 
box radius) from each other in a bounded domain or tend 
to "evaporate" to infinity in an unbounded domain. 
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5 By-product: post-collapse of the 
Smoluchowski-Poisson system 

There is an interesting by-product of the previous study. 
Indeed, the preceding analysis can be used to obtain 
new results concerning the post-collapse dynamics of the 
Smoluchowski-Poisson system (or Kellcr-Scgel model). In 
d = 2 dimensions and for T < T c = GMm/(ik B ), it 
is known that the Smoluchowski-Poisson system forms a 
Dirac peak of mass Mo = (T/T C )M in a finite time t co u 
[9]. For t > t co u the Dirac continues to grow by accretion 
of the surrounding matteiQ. For T = T~ or for t — > +oo 
and any T < T c , the Dirac peak has accreted most of the 
mass. As a result, the system is formed by a Dirac peak of 
mass M(j{t) ~ M surrounded by a dilute halo containing 
the remaining mass e = M — Mo(t). Now, it is possible to 
neglect the self-gravity of the halo and consider that the 
particles of the halo are only subject to the gravity of the 
central peak of mass ~ M. Therefore, the evolution of the 
halo density is governed by a Fokker-Planck equation of 
the form 



k B T. 



-Vp + p- 



GMr 



(165) 



This equation has been studied in [53] in a box. However, it 
was not realized that the corresponding eigenvalue equa- 
tion could be solved analytically in d — 2. Indeed, Eq. 
P65p is equivalent to Eq. (j3~4"]l up to a change of nota- 
tions and we can therefore apply the results of the pre- 
vious sections. It suffices to define D = fcsT/(m^) and 
U = GMmhir. Then, we can use the results of the pre- 
vious sections with now 



(3GMm _ 2T C 
2 ~ ~Y 



(166) 



In an infinite domain, using Eq. (|85p . we find that 
the mass of the Dirac peak saturates to M algebraically 
rapidly as 

^~t- (167) 



1 



In a bounded domain, using Eq. (|149[) . we find that the 
mass of the Dirac peak saturates to M exponentially 
rapidly as 

_ Mg(t) D xl(a)t 



(168) 



where RXi(a) is the first zero of J a -i{x). This result was 
previously found in |24| but the exponential rate (eigen- 
value) was not obtained explicitly (except in the asymp- 
totic limit T -> 0). Note finally that for T = 0, M (t) 
saturates to M in a finite time 



tend — 



R- 



2GM' 



(169) 



4 This post-collapse regime has been studied in [2J] for d > 2. 
For d = 2, it is more difficult to study except in the large 
time limit t — > +oo that we consider here. Note that the virial 
theorem given in [6] is not valid anymore when a Dirac peak is 
formed (in the post-collapse regime). The reason is the same as 
the one given in Sec. 13. II The correct form of the virial theorem 
in that case is given in Appendix [Hi 



corresponding to the deterministic collapse of the outer 
mass annulus initially at r = R (see Eq. (|219[) with k = 
GM/i). 



6 The logarithmic Fokker-Planck equation in 
d dimensions 

In this section, we briefly generalize the previous results to 
the logarithmic Fokker-Planck equation in d dimensions 



dP 
~dt 



p p r 2 



(170) 



The previous results are recovered for d = 2. Taking the 
time derivative of the moment of inertia of the reduced 
particle (|42l) . and using Eq. (|170|) . we obtain after an in- 
tegration by parts 



£1 = -2 



( fc B TVP + PGm^— ) dr. 



(171) 



It can be shown, using the following expressions for P(r, t), 
that the boundary terms at r = and r = +oo vanish. 
Integrating again by parts and introducing the notation 
x(t) = J P(r, t) dr taking into account the possibility that 
the normalization of P(r, t) is not conserved (due to the 
formation of a Dirac peak at r = 0), we obtain 

ti = 2k B Tdx(t)-2Gm 1 m 2 x(t)- (172) 
Finally, this can be rewritten 

l f §=x(*)Mr-T.) ) (173) 
where we have introduced the critical temperature 

M = (174) 
d 

We emphasize that this relation is valid whether P(r, t) is 
spherically symmetric or not. 

We now consider a spherically symmetric evolution. 
Using the notations of Eq. (|35p , we can write the logarith- 
mic Fokker-Planck equation (|170[) in the form 



dP 
~dt 



1 d 



dr 



r d-l 



Dr 



d-i 



dP 
dr 



P 



/3Gmim2 



The time derivative of the normalization is 

J3Gmim2 



f pip 

S d D lim r d ~ x — 
r->o \ dr 



P- 



.(175) 



(176) 



The Fokker-Planck equation (|175[) can be transformed into 
a Schrodinger equation (|49[) with a potential 



V(r) = -§ 



a 2 - 



{d-2f 



(177) 
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where we have defined 

PGm 1 m 2 ~(d-2) 



(178) 



When d > 2, we see that a = at the new critical tem- 
perature 



k B Tl = 



Gm\va 2 
d-2 ' 



(179) 



In the following discussion, it is implicit that T[ = +oo if 
d < 2. The eigenvalue equation (|52l) takes the form 



rV + {d- l)r<f>' + 



A 

J)' 



(d-2) 



= 0, 



(180) 



and it can be solved analytically in terms of Bcsscl func- 
tions. 

Let us first consider the solution 



Mr) 



r^ja^/X/Dr). 



(181) 



Repeating the calculations of Sec. 13.41 we find that the 
solution of the logarithmic Fokkcr-Planck equation (|175l) 
with the initial condition (1551) is 



r + r T 



/ rr \ 



W \2Dt) ■ (182) 



We need to distinguish two cases. For a > (i.e. 
fiGm\m-2 > d — 2 or T < T + '), the behavior of the dis- 
tribution close to r = is 



P+(r,t) 



1 



'0 



1 



2S d r(a + 1) {Dt) a + 1 



e 4Dt 



(183) 



Substituting this equivalent in Eq. (|176p . we find that 



1 



1 



r{a) {Dt)^ 1 



(184) 



In that case, we have an absorbing boundary condition 
at r = and the growth of a Dirac peak. For a < 
(i.e. /3Gmim 2 < d - 2 or T > T*), the behavior of the 
distribution close to r = is 



P+(r,t) 



1 



1 



1 



1 



r f3Gm im2 4 |a| 2S dj T(l + |ffl|) (Dt)\ a \~ 



-e 4Dt . 



(185) 

Substituting this equivalent in Eq. (|176D . we find that 
X+ = 0. In that case, the normalization condition is con- 
served and there is no Dirac peak. 
Let us now consider the solution 



cf>_(r)=r^J_ H (V\/Dr). 



(186) 



Repeating the calculations of Sec. 13.41 we find that the 
solution of the logarithmic Fokker-Planck equation (|175[) 
with initial condition (1551) is 



P-(r,t) 



1 



_L _ra±!f T ( rr \ 

2S d Dt 6 Wt |a| \2Dt) ' 

(187) 



provided that \a\ < 1, according to identity (|70)) . We 
need to distinguish two cases. For < a < 1 (i.e. 
d-2 < f3Gm im2 < d or T* < T < Ti), the behavior 
of the distribution close to r — is 



P-(r,t) 



1 



4" 



1 



r fSGra im2 2S d r{\ - a) (Dty 



-g 4Dt 



(188) 



Substituting this equivalent in Eq. (|176[) . we find that 
X- = 0. In that case, the normalization condition is con- 
served and there is no Dirac peak. For — 1 < a < 
(i.e. d - 4 < /3Gmim 2 < d - 2 or < T < T" with 
ksT'J = Gmirji2/{d — 4)), the behavior of the distribution 
close to r = is 



P- (r, t) 



1 



„d-2 



2\a\ 



1 



1 



2S d r(l-\a\)(Dty-\"\ 



(189) 

Substituting this equivalent in Eq. (|176[) . we find that 
X— > which is not physically possible (the normaliza- 
tion of P(r, t) can decrease if the particle is absorbed at 
r = 0, but it cannot spontaneously increase). Therefore, 
this solution must be rejected. 

These results arc consistent with the van Kampen clas- 
sification of singularities (see Appendix [F| . For T < 
(i.e. a > 1 or /3GtoiTO2 > d), the singularity at r = be- 
haves as an adhesive boundary. In that case, the solution is 
unique and no boundary condition has to be fixed by hand. 
It is given by Eq. (|182[) leading to a Dirac peak (x+ ^ 0). 
For T* < T < Ti (i.e. 0<a<lord-2< /?Gmim 2 < d), 
the singularity at r = behaves as a regular boundary. 
In that case, the boundary condition can be absorbing, 
reflecting or mixed. It has to be fixed by hand. The gen- 
eral solution can be written as a "mixture" (|1 20[) of the 
two solutions (|Tg2j) and (fW)) . The solution (fT52")) leads 
to a Dirac peak (x+ ^ 0) contrary to the solution (|187p 
for which the normalization is conserved (x_ = 0). For 
T > Tl (i.e. a < or /3Gm 1 m 2 < d-2), the singularity 
at r = behaves as a natural repulsive boundary. In that 
case, the solution is unique and no boundary condition 
has to be fixed by hand. It is given by Eq. (|182D which 
does not lead to a Dirac peak (x+ = 0). 



7 Conclusion 

In this paper, we have analytically studied the evolution of 
N = 2 Brownian particles in gravitational interaction in 
a space of d = 2 dimensions. Up to a change of notations, 
this is equivalent to the simplified motion of two biologi- 
cal entities interacting via chemotaxis (in which case the 
dimension d = 2 is physically relevant). Of course, the 
consideration of only N — 2 particles is an extreme limit 
but the problem is already involved and shows that the 
dynamics is complex since the particles can coalesce to 
form Dirac peaks. The same phenomenon (collapse and 
Dirac peaks) occurs for a larger number of particles and 
has been investigated analytically in the mean field limit 
N — > +oo [819110] . It shares some analogies with the Bose- 
Einstcin condensation [25] . The case of a finite number of 
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particles will be investigated numerically in a forthcoming 
paper [TT] . 

Finally, we would like to point out some analogies 
with the transport of passive particles by a stochastic 
turbulent flow characterized by scale invariant structure 
functions (Kraichnan model) [26 27 28 or, more gener- 
ally, with correlated Brownian motions with scale invari- 
ant correlations |29j . In particular, implosivc collapse of 
trajectories has been found by Gawedzki & Vergassola [2H1 
for strongly compressible flows. These authors determined 
the statistics of inter-trajectory distances and observed a 
lack of normalization when the diffusivity tends to zero. 
Like in our problem, a defect of probability concentrates 
at r — in a (^-function term carrying the missing prob- 
ability. These authors also studied the long time behav- 
ior of the averaged powers of the distance between the 
Lagrangian trajectories and obtained power law behav- 
iors. In these hydrodynamical problems, the choice of the 
boundary condition at r = when viscosity and diffusiv- 
ity go to zero is crucial and different regimes have been 
investigated. In the case of weak compressibility, the sin- 
gularity at r — acts as a repulsive entrance boundary 
and the particles never collide. For strong compressibility 
and smooth flows, the point r — behaves as a natural 
attractive boundary for which particles approach r = in 
an infinite mean time. For strong compressibilities, r = 
works as an adhesive boundary for which particles col- 
lide in finite time with a vanishing relative velocity. In 
that case, they remain at that point indefinitely (coales- 
cence). Since the probability to find at t > the particle 
at r — is finite and increases in time, no stationary state 
is reached and P(r, t) develops a Dirac delta function at 
r = with a time increasing coefficient. For intermedi- 
ate compressibilities, r = works as a regular boundary 
because particles hit one another in finite time but with 
non-zero relative velocity. In that case, both attractive and 
repulsive solutions are possible and it is necessary to fix 
by hand a boundary condition at r = 0. These different 
regimes have been discussed in these terms by Gabriclli & 
Cecconi [29] in relation to boundary (or singularity) clas- 
sification introduced by Van Kampen [3U] (see also Feller 
|31j). A similar phenomenology is obtained in the framer- 
work of our Brownian model. 

Acknowledgment: We thank the referee for mentioning 
the connection of our study with the transport of a passive 
particle by a stochastic turbulent flow (Refs. |26I27I28I29) ) 
and indicating to us the van Kampen classification. This 
led to a more detailed analysis of our model. 



A Moment of inertia and critical 
temperatures 

It is shown in (see also Appendix QI that, in d — 2, 
the total moment of inertia of self-gravitating Brownian 
particles 



If.nt — 



{™<*r 2 a ), 



(190) 



satisfies the virial theoren0 
1 Ahot 



Nk B (T-T c )-PV, 



(191) 



where P is the total pressure at the boundary of the do- 
main. The critical temperature appearing in this relation 
is 

= GE a ^m a m, = G_( M2 _j : ,\ 
AN 4N y / 

For equal mass particles 



k B T c = (N-l) 



Gm 2 



(193) 



In fact, it may be more relevant to measure the posi- 
tions of the particles relative to the center of mass 



R 



M 



and define the moment of inertia by 

/ = ^](?7i Q (r Q -R) 2 ) 



(194) 



(195) 



Indeed, if all the particles collapse in a single point, this 
point will be the center of mass R. Therefore, / will be 
zero while I tot is non zero. We clearly have the relation 



M(R 2 



(196) 



We also check that, for N = 2, I represents the moment 
of inertia of the reduced particle. The center of mass has 
a pure Brownian motion 



with a diffusion coefficient 



k B T 



It satisfies a relation of the form 



(197) 



(198) 



(199) 



where P e ff is an effective pressure on the boundary of 
the box due to the center of mass. Therefore, the virial 
theorem expressed in terms of I is 

I^f" = Nk B(T-T c )-k B T-APV, (200) 

where AP = P — P e ff- Since the motion of the center 
of mass is completely decoupled, it is as if we had only 



According to the discussion of Sec. 12.41 we now know that 
this relation ceases to be exact when the particles form Dirac 
peaks. However, we shall not address this problem here. 
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N — 1 particles in the system. Therefore, it makes sense 
to rewrite the foregoing equation in the form 



y^ = (N-l)k B (T-T.)-APV, 



where 



N 



-knTr. 



N - 1 

For equal mass particles, we have 



(201) 
(202) 

(203) 



Furthermore, if we consider an infinite domain and mea- 
sure the displacement of the particles relative to the cen- 
ter of mass, using Eq. (|201|) . we find that ((r — R) 2 ) = 
AD(T)t+ ((r — R) 2 )o with an effective diffusion coefficient 



. . N-l k B T f T* 
D(T) = ^— 1 - 



(204) 



For N = 2 particles, noting that /i = to/2 and ((r— R) 2 ) = 
(((n-R) 2 ) + ((r 2 -R) 2 ))/2 = (r 2 )/4, where (r 2 ) denotes 
the mean square displacement of the reduced particle, we 



find that Eq. (j204|) is consistent with Eq. (|30|). 

Coming back to the general expression (|201[) . we con- 
clude that the system is expected to collapse and form 
Dirac peak(s) for T < and evaporate (in an infinite do- 
main) or tend to an equilibrium state (in a finite domain) 
for T > T*. The temperature T* is precisely the collapse 
temperature that was obtained in |12| directly from the 
study of the partition function in a bounded domain (the 
partition function diverges for T < T*). We now under- 
stand better the relationship between the two tempera- 
tures T c and T*. The temperature T* is associated to the 
collapse of the system. It arises in the virial theorem after 
the contribution of the center of mass has been removed. 
However, the critical temperature appearing in the exact 
equation of state is T c . Indeed, at statistical equilibrium, 
wc have the relation 1121: 



PV = Nk B {T-T c ), 



(205) 



which also results from Eq. (|19ip . At T = T*, the particles 
form a Dirac containing all the mass and the equation of 
state (|205p reduces to 



PV = fc B T» 



(206) 



where we have used Eq. (|202j) . This is the equation of state 
of a single Brownian particle (the Dirac) at temperature 
T = T*. It is of the usual form PV = Nk B T with N = 1. 



B Check of consistency for T ^ 

In this Appendix, we check that the relation (|92p is con- 
sistent with the virial theorem (|43[) . The virial theorem 
(l43l) can be rewritten 



i = ADfi X (t){l-a). 



(207) 



Integrating this relation, we get 



I(t) = AD(1 - a)(i / X (r)dT + I(0) 



(208) 



with 1(0) = A/r 2 in our case. Using Eq. (|78j) . we have 



X(r)dr = t- / r \ ^ dr. 
10 Jo r(a) 

Setting x — r 2 / (4Dr) , this can be rewritten 

/"* / w 1 T l f + °° r, , s dx 

x ( T )dT = t-—J- l 2 r a (x)—. 

Jo r(a) 4D J la x 2 



Let us consider the integral 



r a (x) 



dx, 



K(s) -- 
where we recall that 

r a (x) = I t a - L e- T dt. 

J X 

Integrating by parts, we get 



K(s) 



r a (s) 



x a ~ 2 e- x dx. 



Integrating by parts again, we obtain 



K(s) =(- + 



s 1 — a 

Wc now have 
I(t) = iirl+4D(l-a)n t 



r a (s) - 



1 - a 



1 rl K (^_\ 



(209) 



(210) 



(211) 



(212) 



(213) 



(214) 



(215) 



r{a)AD \ADt J 
Substituting Eq. (12141) in ([215]) . we recover Eq. ([92]). 

C The case T = 

The case T = can be treated specifically. In that case, 
the reduced particle has a deterministic motion given by 
the equation 

dr r , 

*=-*;*' (216) 



where 



k = 



GTO1TO2 

— 1^' 



(217) 



This equation can be integrated at once. If a denotes the 
initial position of the particle, its position at time t is 



r 2 = a 2 - 2kt. 



(218) 



This study can be easily generalized in d dimensions. 
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Therefore, the particle reaches the origin at a time 



(219) 



Equivalently, the particle that reaches the origin at time 
t was located initially at 



a(t) 



<2kt. 



(220) 



Let Po(a) be the initial probability density to find the 
particle in a. The conservation of the probability density 
imposes 

P(r, t)2irrdr = P (a)2irada. (221) 

Now, according to the equation of motion (|218p . we have 
rdr = ada so that P(r,t) = Pq(ci). Therefore, the proba- 
bility density to find the particle in r at time t is 



P(r,t)=P (Vr 2 + 2kt) 



(222) 



We can check by direct substitution that this is indeed the 
solution of the Fokkcr-Planck equation (|34|) at T = 0: 



(223) 
= at 

(224) 



dP _ kdP 
dt r dr 

The probability that the particle has not reached r 
time t is 



xit) 



P(r,t)2irr dr. 



Using the distribution (|222[) and performing the change of 
variables (|218p . we get 



x(t) = 



/-foo py2kt 
_P Q (a)2irada=l- P Q {a)2ira 
/2fct JO 



da. 



(225) 



Therefore, the probability that the particle has formed a 
Dirac peak at time t is 



pV2kt 

XD (t) = / P (a)27ra 
Jo 



da. 



(226) 



This corresponds to probability to find initially the par- 
ticle in the disk of radius y/2kt. For consistency, let us 
derive this result in a different manner. Starting from the 
relation 



XD =2irkP(0,t), 
and using Eq. (|222|) . we get 

X D = 27TkP Q (VWt) . 

Integrating this relation, we find that 



(227) 



(228) 



XD{t) = 2irk / P 
Jo 



'2kr\ dr 



/2k I 



Po(a)2na da. 

(229) 



This returns Eq. (|226p as it should. Finally, the moment 
of inertia is 



r+oo 

I(t) = / P(r,t)iir 2 2nrdr. 
Jo 



(230) 



Substituting the probability density (|222[) in Eq. (|230j) and 
performing the change of variables (|218p , we obtain 



I(t) 



+00 



Po(a)/i(a 2 — 2kt)2-Kada. 



2kl 



This can be written equivalently 



/+oo 
_ P {a)/j,a 2 2na da - 2knx(t)t- 
,l2kt 



(231) 



(232) 



In Appendix [D] we check that this relation is consistent 
with the virial theorem (|43D . 

Let us specifically apply these results to an initial dis- 
tribution of the form (|60p . Since the motion is determin- 
istic, the particle initially located at tq will be located at 
r(t) = — 2kt at time t. It will reach the origin in 
a finite time t c = rp/(2fc) (at that time, the two original 
particles stick together and remain tightly bound) . There- 
fore, the probability that the particle has formed a Dirac 
at time t is a Heaviside function: xd = if t < t c and 
Xd = 1 if t > t c . The moment of inertia is I(t) = ^r(t) 2 
if t < t c and 1(0) = if t > t c . 



D Check of consistency for T = 

For T = 0, the virial theorem (|4"3")) becomes 
Integrating this relation, we get 



I(t) = 1(0) - 2kn / X (r) dr. 

Jo 



Integrating by parts, we have 



x(t) dr = tx(t) - / x(t)t dr. 
o Jo 



According to Eqs. ([39]) and (|222p . we have 
X(r) = -27rfcP (v / 2fcI). 

Therefore 



(233) 



(234) 



(235) 



(236) 



= 1(0) - 2knx(t)t - Airk^fi / P (V2kT) T dT. (237) 

Jo 

Setting a = \f2kr, we obtain 

I(t) = 1(0) - 2kvx(t)t - 2717* / Po(a)a 3 da. (238) 

Jo 



Since 



1(0) = 2-Kfi / P a (a)a 6 da, 



(239) 



finally recover Eq. (I232D . 
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E Asymptotic behaviors 



F The van Kampen classification 



In this Appendix, wc determine the asymptotic behaviors 
of the function 



i 



[xr{x)]v*' 



(240) 



appearing in Eq. (j8"T)) . We first determine the behavior of 
the Gamma function for small x. Expanding the identity 

r(i - x)r(x) 



sin(7ra;) ' 
for x — > 0, and using 

r'(l) = -7 ~ -0.57721..., 
where 7 is the Euler constant, we obtain 

xT(x) ~ 1 — <yx, (x — > 0). 
From the previous result, we deduce that 

\n{[xr(x)] 1/x } ~ ~ln(l-~/x) -> 

Therefore, for x —> 0, we find that 
1 



-> e> 



1.78107. 



[xr(x)Y/ x 

On the other hand, using the equivalent 

r(x) ~ %/27re~ :E .T a:: "5 j (3; _> +0O ) ; 

we find, for x — > +00, that 

1 e 
[iffa;)] 1 / 1 ~ sc' 

Considering now the function defined by Eq. 
using the equivalent 

r(aO~i, (x->o), 

x 

we find, for a — > 0, that 

- aP (l) - 0.219384a. 

To determine its asymptotic behavior for a — > 
first note that r a (x) = r(a) — 7(0,2;) where 



7(0,0;) = / e~H° 
Jo 



For a 



- x dt. 

-00, we have the asymptotic expansion 

7(0,2;) 



^ (a + n)n\ 

n— 



(241) 

(242) 
(243) 

(244) 
(245) 

(246) 

(247) 
and 

(248) 

(249) 
+00, we 

(250) 
(251) 



This implies 7(0, 1) ~ l/(oe). On the other hand, 

r{a) - V2^ e - a a a -^. (252) 

Combining these results and recalling that Xd(1) = 1 — 
7(0, l)/r(a), we obtain 



1 



D a-t 



(253) 



In this Appendix, we apply to our system the boundary 
classification introduced by van Kampen |30j . For a sum- 
mary, we refer to Gabrielli and Cecconi [53]. In order to 
avoid repetitions, we shall use their notations and we refer 
the reader to their paper for more details. 

For the sake of generality, we consider the logarithmic 
Fokker-Planck equation in d dimensions 



dP 
~&t 



1 d 



Dr 



d-i 



dP 

dr 



P 



f3Gmim,2 



.(254) 



To apply the van Kampen classification, we must first 
transform the Fokker-Planck equation (|254[) into a one 
dimensional Fokker-Planck equation. To that purpose, we 
set r = \[Dx and f(x,t) = VDP(r,t)S d r d - 1 . This trans- 
forms Eq. (|254j) into 



dl = d_(d£ 
dt dx \ dx 



J. + f3Gm 1 m 2 , 



(255) 



with the normalization condition J Q f(x, t) dx = 1. This 
is a one dimensional Fokker-Planck equation of the form 

with D(x) = 2 and K(x) = (d — 1 — ^Gm\m2)/x. Van 
Kampen's classification for a singularity x = is based on 
the analysis of the behavior for e — > of the integrals 



U = 



dx e^ x) 



, r x e -<t>{*') 

him- 



l. 



dx 



-<Kx) 



where 



(x) = -2 / da; 



D(x') ' 



(257) 



(258) 



(259) 



(260) 



and with xq > 0. For the logarithmic Fokker-Planck equa- 
tion (|255p . we have D(x) = 2 and 



(261) 



4>(x) = —(d — 1 — /3GmiTO2) In — 

\a;o 



for a —} +00. 



Considering the limit e — > 0, it is easy to see that: (i) L\ < 
+00 iff f3Gm\m2 > d — 2, (ii) L2 < +00 iff /3Gmim2 > 
d—2, (iii) L3 < +00 iff f3Gmini2 < d. Therefore, according 
to van Kampen's classification, we need to consider three 
cases (it is useful to introduce the critical temperatures 
= Gm\m,2i 'd and T + ' = Gmim2/(d— 2)): 
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(i) if (3Gmim 2 < d - 2, i.e. T > T*, the singularity 
x = behaves as a natural repulsive boundary (L\ — > 
+oo). The particle run away from the singularity never 
touching it. The solution is unique. 

(ii) if d - 1 < f3Gm im2 < d, i.e. < T < T^, 
the singularity x = behaves regular boundary 
(L\, L2, L3 < +00). In that case, an absorbing or reflecting 
boundary condition has to be fixed by hand to determine 
the solution of the equation. 

(hi) if [3Gmim,2 > d, i.e. T < T„, the singularity x = 
behaves as an attractive adhesive boundary (L\, L 2 < +00 
and L3 — > +00). In that case, f(x, t) develops a Dirac peak 
at x = with a time increasing coefficient. The solution 
is unique. 

For d < 2, we only have a transition at temperature 
= Gmiiri2/d. For d > 2, we have two transitions at 
temperatures = Gm\m,2l d and = Gm\m2j{d — 2). 



G Temporal correlation functions and front 
structure of the logarithmic Fokker-Planck 
equation 



We consider the logarithmic Fokker-Planck equation (|254l) 
in a space of dimension d. We assume that the do- 
main is unbounded. In order to have an equilibrium state 
P e (r), the potential must be regularized at short distances. 
Therefore, we assume that the potential has a logarithmic 
behavior for sufficiently large r and that it tends to a fi- 
nite constant for r —> 0. In that case, P e oc r -P G ™i™2 
for r — > +00, and there exists an equilibrium state iff 
fcflT < fcflT, = Gm\m2l d. To determine the temporal 
correlation functions, we use the theory of Marksteiner 
et al. [32] (see also 133134135136137138] ). As shown in Ap- 
pendix [F] the Fokker-Planck equation (|254p can be trans- 
formed into a one dimensional Fokker-Planck equation of 
the form 



dt dx \ dx dx 



(262) 



with a potential behaving like <P(x) ~ a In 2 for x — > +00, 
where a = fiGm\m<i — (d— 1). We introduce the temporal 
correlation functions C(t) = (A(0)A(t)) — (A) 2 and refer 
to |32I33I34I35I36I37I38] for the details of the calculations 
(we use here the notations of [37]). If A(t) = x(t) n , then 



C(t) - t 



t = -ri- 



ot - 1 



(263) 



In the present case, the exponent can be written £ = 
— n + {fiGm\m2 — d)/2. Equation (|263[) is valid provided 
that £ > 0, i.e. fcgT < GraiW2/(2n + d), which corre- 
sponds to the condition of existence of the moment (x 2n ) 
at equilibrium. If A(t) ~ (lnx(t)) 1/s , then 



C(t) 



(Int) 



2/S 



t~ 



(264) 



We note, finally, that the relaxation of the tail of the 
distribution function P(r, t) that is solution of the loga- 
rithmic Fokker-Planck equation (|254l) can be studied with 



the approach developed by Chavanis & Lemou [35]. In 
particular, the function u(r,t) = P(r,t)/P e (r) has a front 
structure and the position of the front evolves in time like 
r f (t) - \J1Dat. 

H Virial theorem in the post-collapse regime 
of the Smoluchowski-Poisson system 

A virial theorem associated with the Smoluchowski- 
Poisson system in d = 2 dimensions has been derived in 
[5]. However, the derived equation is not valid in the post- 
collapse dynamics when a Dirac peak forms at r = and 
grows. This happens for T < T c = GMm/ '(4fee) when 
t > t co u [B]. The reason is the same as the one given in 
Sec. 13.11 In this Appendix, we provide the proper form of 
the virial theorem that is valid both in the pre and post 
collapse regimes. 

The total density profile of the self-gravitating Brow- 
nian gas can be written 



p(r,t) = M D (t)S(r) + p(r,t), 



(265) 



where the first term takes into account the possible for- 
mation of a Dirac peak at r = and the second term is 
the (regular) density profile excluding the Dirac. The to- 
tal mass is M = M D (t) + M (t) where M D (t) is the mass 
contained in the Dirac peak and M(t) = J p(r, t) dr is the 
mass outside the Dirac. The Smoluchowski-Poisson sys- 
tem accounting for the presence of a Dirac peak can be 
written 



dp 
^dt 



k B T^ GM D (t) 



A<P = S d Gp. 



pW<P 



(266) 



(267) 



Using Mb = —AI and integrating Eq. (|266|) on the infinite 
space, we find that the mass accumulated in the Dirac 
peak by unit of time is 



dM 



D 



dt 



SdG 



M D (t)p(0,t). 



(268) 



Equations (|266l) - (|268[) form a closed system describing 
the evolution of the system in the pre and post collapse 
regimes. These equations have been studied in fM\ . 

We now specialize on the 2D Smoluchowski-Poisson 
system. Taking the time derivative of the moment of in- 
ertia 



pdr, 



(269) 



and using Eq. (|266|) . we obtain after integrations by parts 

dJ = 4kBT M( Q _ 2GM{t)M D (t) + 2W U , (270) 
at m 

where Wu = — J pr ■ V<P dr is the virial of the gravita- 
tional force. In d = 2 dimensions, it is equal to Wu = 
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-GM(t) 2 /2 (see [5J and Appendix | 
tain the virial theorem 



Therefore, we ob- 



k B T GM D (t) GM(t) 



(271) 



that is valid in all the regimes of the dynamics. For T > T c 
or in the pre-collapse regime t < t co ii for T < T c , there 
is no Dirac peak at r = 0. In that case, Mjj(t) = 0, 
M(t) = M, and the virial theorem (|27ip reduces to 



yi t =Nk B (T-T c 



(272) 



where k B T c = GMm/A. This returns the result of [B]. 
Let us now consider the post collapse regime t > t co u for 
T < T c . When t —> +oo, almost all the mass is in the 
Dirac so that M D (t) ~ M and M(t) = M - M D (t) ~ 
0. In that case, the Smoluchowski equation (|266[) can be 
approximated by 



dp 



kuT. 



m 



GAI 



and Eq. (|268[) becomes 

dM D S d G 



dt 



i 



Mp(0,t). 



(273) 



(274) 



This is equivalent to Eqs. (|34D and (|4"U)) studied in this 
paper, with a simple change of notations discussed in Sec. 
[5j In that case, the virial theorem (|27ip becomes 



(275) 



where k B T* = GMm/2 (i.e. T* = 2T C ). This is equivalent 
to the virial theorem P5|) . 



I Virial theorem for power-law interactions 

In this Appendix, we provide the general form of the virial 
theorem for Brownian particles with power law interac- 
tions in d dimensions. We only give the final expressions, 
and refer to [5J for more details on their derivation. As 
explained in Sec. 12.41 the following expressions are valid 
as long as there are no Dirac peaks. 

Let us consider N Brownian particles with individual 
mass m a in a space of dimension d. We assume that the 
particles are subject to an external harmonic potential 
V(r) = \u 2 r 2 and that they interact through an alge- 
braic potential = — d+7 „ 2 if 7 7^ 2 — d and 
through a logarithmic potential u(£) = Gln£ if 7 = 2 — d, 
both corresponding to a force —u'(£) = — G/£ +7_1 . The 
gravitational potential is recovered for 7 = 0. The case 
studied in the present paper is very particular because it 
corresponds to a logarithmic (7 = 2 — d) and a Newtonian 
(7 = 0) interaction. The stochastic equations of motion of 
the particles are 



ld+7 



£c? + y/2D2B? (t), 
(276) 



where B Q (i) is a white noise. Here, the Greek letters refer 
to the particles and the Latin letters to the coordinates 
of space. The diffusion coefficient is given by the Einstein 
formula D a = ^k B T/m a . The moment of inertia tensor is 
defined by 

J«=5^m a i?a:?. (277) 



We introduce the kinetic energy tensor 



and the potential energy tensor 



Wij = G m a 



r /3 



ld+7 



1 (s? -af 



where the second equality results from simple algebraic 
manipulations obtained by interchanging the dummy vari- 
ables a and /? and summing the resulting expressions. The 
tensor virial theorem associated with the stochastic equa- 
tions (|2"75]t is 



1 



1 



2 13 2 s 13 

1 



(P ik Xj + PjkXi) dS k , 



(280) 



where the last term takes into account pressure forces at 
the boundary of the system. For Brownian particles, it is 
implicitly assumed that the quantities appearing in Eq. 
(|280[) are averaged over the noise and over statistical real- 
izations, while for Hamiltonian systems (£ = D a = 0), Eq. 
(|280p is exact without averages. The scalar virial theorem 
is obtained by contracting the indices leading to 



P ik XidS kl (281) 



where 



Tfl a X n 



K = \ £ m 



1 



(282) 



are the moment of inertia and the kinetic energy. On the 
other hand, Wu is the virial which takes the form 



For 7 7^ 2 — d, we find that 

W u = (d + 7 - 2)W, 
where W is the potential energy 

G \—\ m a mp 



W = - 



£ 



2(d + 7 -2)^|r /J -r a |*W-a- 



(283) 



(284) 



(285) 
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In that case, the scalar virial theorem reads 



1+^1 + uj%I = 2K + (d + 7 - 2)W - <f> P ik XidS k . 



(286) 



For Hamiltonian systems (D = £ = 0), the scalar virial 
theorem in an unbounded domain (P = 0) reduces to 
[10131] : 



-I + wgl = 2K + (d + 7 - 2)W. 



(287) 



Since the total energy E = K + W + is conserved, 

we obtain 



-I + 2ujp = 2E + (d + 7 - 4)W. (288) 



For the index 7 = 4 — d, we get 



I + AljZI = 4P. 



(289) 



When Wq = 0, we obtain I = AE which yields after in- 
tegration I = 2Et 2 + C\t + C 2 . For E > 0, I +00 
indicating that the system evaporates. For E < 0, I goes 
to zero in a finite time, indicating that the system forms 
a Dirac peak in a finite time. When ui 2 > 0, the moment 
of inertia I oscillates with pulsation 2luq around the value 
E/uj 2 . For uj 2 = —^0 < 0' corresponding to a repulsive 
harmonic potential (or a rotation) , the moment of inertia 
increases exponentially rapidly as e 217 °*. For the gravita- 
tional interaction (7 = 0), the relation (|289p is valid in a 
space of dimension d — 4 [5]. For d = 1, this relation is 
valid for the potential u = — This is related to the 

Calogero-Sutherland model [42143] . 

For 7 = 2 — d, corresponding to a logarithmic potential 
in d dimensions, we have the simple exact result 



--G ^ m a mp. 



(290) 



It is interesting to note that this expression only depends 
on the mass of the particles and not on their position. For 
equal mass particles, 



Since 



W H = ~GN(N -l)m 2 . 



}^ m a mp = M 2 - m 2 a , 



(291) 



(292) 



we see that the first term is of order N 2 m 2 and the second 
of order Nm 2 (where m is a typical mass). Therefore, in 
the mean- field limit N — > +00, we obtain 



IF 



m.f. _ 



GM 2 



2 ' 



(293) 



whatever the number of species in the system. 



At equilibrium, the scalar virial theorem (|281[) reduces 



to 



2K + Wu -up= I PikXidSk 



(294) 



For Hamiltonian systems, this relation is valid for a steady 
state after time averages, or averages over statistical re- 
alizations, have been made. If the system is at statisti- 
cal equilibrium, then K = ^NksT and Py = pSij with 
p = J2 S PskBT/m s , where p s refers to the density of the 
different species. Introducing the notation P = § pr-dS 
[6], we get 



dNk B T + W H 



dPV. 



(295) 



For an ideal gas without interaction (Wu — 0), we recover 
the perfect gas law PV + uj 2 I/d = NksT in the presence 
of a harmonic potential (when lu 2 > 0, we get in an un- 
bounded domain / = dNksT/uiQ, and when u>q = 0, we 
get PV = NksT). Alternatively, for a gas with logarith- 
mic interactions (7 = 2 — d), using Eq. (|290[) . we obtain 
the exact equation of state 



PV+^-=Nk B (T-T c ), 
with the exact critical temperature 

kBTc = 2dN • 

For equal mass particles, we get 

k B T c = (N-l)- 



2d 



In the mean-field limit 



k B Tr l = 



GM 2 
2dN 



If wo =0, the equation of state (|296p reduces to 
PV = Nk B (T-T c ). 



(296) 



(297) 



(298) 



(299) 



(300) 



When u 2 > 0, according to Eq. (|296[) . an equilibrium state 
can possibly exist only for T >T C (since P > and I > 0). 
For T = T c , we have P = if uj'I = 0, P = I = if uj 2 > 
and / = dPV/f2 2 if ujg = — J?g < 0. In an unbounded 
domain (P = 0), we get 



Nk B (T — T c ). 



(301) 



If ujq = 0, an equilibrium state can possibly exist only 
for T = T c . If uiq = — i?Q < 0, an equilibrium state can 
possibly exist only for T < T c . For T = T c and ojq ^ 0, we 
must have 7 = 0. 

We now consider the strong friction limit £ — > +00 
where inertial effects are negligible. In that limit, the 
velocities thermalize on a timescale of order l/£. In 



that case, K; 



^Nk B TSij and P^- 



pdij with p 
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Ps^bT /m s even if the system has not yet reached a 
state of mechanical equilibrium [6] . From Eq. (|280[) , we ob- 
tain the overdamped virial theorem for a self-gravitating 
Brownian gas 



-Zhj + wlhj = Nk B T8 l3 + Wij 



- J p(xidSj + XjdSi). 

(302) 

We can obtain this result in a different manner. In the 
strong friction limit £ — > +oo, the inertial term in Eq. 
(|276p can be neglected so that the stochastic equations of 
motion reduce to 



Gm^-xf) <4 



VP 



(303) 

where D' a = k B Tfi a is the diffusion coefficient is physical 
space and /i a = l/(£m a ) the mobility. The overdamped 
virial theorem (|302j) can be directly obtained from these 
stochastic equations [B]. The scalar virial theorem reads 



^1 + cup = dNk B T + Wu 



dPV. 



For 7 = 2 — d, using Eq. (|290[) . wc obtain 

hi + wgj = dNk B (T - T c ) - dPV, 



(304) 



(305) 



with the exact critical temperature (|297p . In an infinite 
domain (P = 0), this relation reduces to 



-M + Lj 2 Q I = dNk B (T-T c ). 



(306) 



This is a closed equation that can be solved analytically. 
For wq 7^ 0, the solution is 



I(t) 



1(0) 



dNk f 



(T - T c ) 



dNk B , 

2^(T-T C ). 



Let us introduce the new critical temperature 

.2^(0) 



kuT us — kftT r 



dN 



(307) 



(308) 



depending on the initial value of moment of inertia (note 
that T u > T c ). This is the value at which the term in 
bracket in Eq. (|307[) vanishes. If T > T c , the system tends 
to an equilibrium state corresponding to I eq = dN Ji B (T — 

T c ), sec Eq. poip . More precisely, for T > T u , the moment 
of inertia increases, for T c < T < T u it decreases and 
for T = T u it remains constant. If T = T Cl wc find that 
I(t) —> for t —> +oo implying a collapse in infinite time. 
If T < T r , the moment of inertia vanishes at a time 



In 



dNk B (T c — T) 



(309) 



implying the finite time collapse of the system (we recall 
that I(t) = M(r 2 ){t)). For uj = 0, wc obtain 



I(t) 



2dNk B 



(T-T c )t + I{0). 



(310) 
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Fig. 12. Time evolution of the moment of inertia for an over- 
damped Brownian system with logarithmic interactions and 
attractif harmonic potential. 
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Fig. 13. Time evolution of the moment of inertia for an over- 
damped Brownian system with logarithmic interactions. 



If T > T c , the system evaporates. If T = T c , wc find that 
I(t) = 1(0) for all times. If T < T c , the moment of inertia 
vanishes at a time 



Lend 



^(0) 



2dNk B (T c — T) ' 



(311) 



implying the finite time collapse of the system. If ujq = 
— f?Q < (repulsive harmonic potential or rotation), the 
picture is different. For T > (note that now T u < T c ), 
the system evaporates. For T < T u , which is possible iff 
1(0) < dNk B T c J (— u! 2 ), the moment of inertia vanishes 
at a time t en d, given by Eq. (|309[) . implying finite time 
collapse. Finally, for T = T u , the moment of inertia is 
conserved. Some representative curves of these different 
evolutions are given in Figs. [T2l [TBI and IT4l 

Wc now give the proper form of virial theorem corre- 
sponding to the generalized Smoluchowski equation [B]: 



3e = v- 

dt 



- (Vp + + poj%r) 



(312) 



where the pressure p = p(r, t) is given by an arbitrary 
barotropic equation of state p = p(p). We assume that 
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Fig. 14. Time evolution of the moment of inertia for an over- 
damped Brownian system with logarithmic interactions and 
repulsive harmonic potential. 



the particles are subject to an external harmonic potential 
V(r) = ^0J 2 r 2 and that they interact through an algebraic 



potential 



#(r,t) = - , ° - / , ^§ , dr', (313) 
v ; rf + 7-2 7 | r _r'|<i+7-2 v > 

if 7 7^ 2 — d, or a logarithmic potential 

<P(r,i)=G / p(r',t) In |r - r'| rfr', (314) 

if 7 = 2 — d. The mean force of interaction acting on a 
particle in r is 

F = — W = ~G J P(r\t) lr r _-^ dr'. (315) 

For simplicity, we assume that the particles have the same 
mass to. The Lyapunov functional associated with the 
Smoluchowski equation (|312[) is the free energy 

F = J p d P ' dr +\ J P®dr + J pVdr, (316) 

and it satisfies an iJ-theorem, i.e. F < 0. The Smolu- 
chowski equation (|312j) with an isothermal equation of 
state p(r,t) = p(r,t)kBT/m is the mean field Fokker- 
Planck equation associated with the overdamped stochas- 
tic process (|303[) . The potential energy tensor is defined 
by 

d<P 



Wi, 



px t —— dr, 

OX 4 



while the virial is 



W H = - pr- \7<Pdr. 



(317) 



(318) 



Substituting Eq. (|315[) in Eq. (|3 1 7[) and using the usual 
symmetrization procedure, it can be rewritten 



Wij = -rG / p(r)p(r') 



-^i)i x i - x 'j) 



|r — r'\ d+ ^ 



dvdv'. (319) 



Contracting the indices, we get 

^. = 4 G /ji^p&**' < 320 > 

For 7 7^ 2 — d, we obtain 

W u = (d + 7 - 2)W, (321) 

where W = \ f p<Pdr is the mean field potential energy. 
For 7 = 2 — d, we obtain 



Introducing the moment of inertia tensor 



pXiXj dr, 



(322) 



(323) 



we find that the tensor virial theorem associated with the 
generalized Smoluchowski equation (|312[) is given by 



1 ... 
2 



+ ui%Iij = 5ij J pdr + Wij - ^j> p(xidSj + XjdSi). 

(324) 

The scalar virial theorem, obtained by contracting the in- 
dices, takes the form 

f J + top = d j P dr + W u - dPV, (325) 



1 ^ 

2 



where 



1=1 pr 2 dr, 



(326) 



is the moment of inertia. The equilibrium scalar virial the- 
orem is 



wg J = d / pdr + Wii-dPV. 



(327) 



For an isothermal equation of state p = pksT /m, we re- 
cover Eq. (|304[) where W^j is now given by Eq. (|322p . For 
a logarithmic potential (7 = 2 — d), we recover Eq. (|305D 
where T c is given by 



kf)T r 



GNm 2 
2d 



(328) 



At equilibrium, we recover Eq. (|295l) . 

Finally, we give the proper form of virial theorem for 
the damped barotropic Eulcr equations 



dp 

at 



v • (pu) = 0, 



(329) 



o -I 

— + (u- V)u = --Vp- V<£-£u-cjgr, (330) 

under the same conditions as before. The tensor virial the- 
orem is given by 



2 ItJ + 2' 



Uj + + Wglij = I pu, 11 , dr + 5ij j pdr 



~-f p(xidSj + XjdSi), (331) 
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and the scalar virial theorem by 

^1 + + UqI = j pv? dr + d J pdr + W u - dPV. 

(332) 

At equilibrium, we obtain Eq. (|327|l . The virial theorem 
for the barotropic Eulcr equations is recovered by taking 
£ = |4QI41j . The Lyapunov functional associated with 
the damped Euler equations (|329[) - (|330p is the free energy 

F = J p J P ^-dp'dr+^ J p$dr+ J P Vdr+ J p\dr, 

(333) 

and it satisfies an iJ-theorem, i.e. F < if £ ^ 0. For 
the Euler equations (£ = 0), the energy functional (|333D 
is conserved F = 0. 
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